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This paper introduces a spatio-temporal resonator model and an inference method for detection 
and estimation of nearly periodic temporal phenomena in spatio-temporal data. The model is de- 
rived as a spatial extension of a stochastic harmonic resonator model, which can be formulated in 
terms of a stochastic differential equation (SDE). The spatial structure is included by introducing 
linear operators, which affect both the oscillations and damping, and by choosing the appropriate 
spatial covariance structure of the driving time- white noise process. With the choice of the linear 
operators as partial differential operators, the resonator model becomes a stochastic partial differen- 
tial equation (SPDE), which is compatible with infinite-dimensional Kalman filtering. The resulting 
infinite-dimensional Kalman filtering problem allows for a computationally efficient solution as the 
computational cost scales linearly with measurements in the temporal dimension. This framework 
is applied to weather prediction and to physiological noise elimination in fMRI brain data. 

PACS numbers: 82.40.Bj, 89.75.Kd 



I. INTRODUCTION 

Oscillations stem from repetitive variation, typically in 
time, of some measure around a point or an equilibrium. 
This type of phenomenon is commonly encountered in 
natural systems, as well as in physical, biological, and 
chemical models [1-3]. This paper proposes a computa- 
tionally effective evolution-type stochastic partial differ- 
ential equation model and an inference method, which 
together provide a novel and efficient means of detecting 
and modeling latent oscillatory structures in space-time, 
such as physiological noise in fMRI brain data [4, 5] or 
temperature variation in climate models. 

The proposed model can be thought of as an exten- 
sion of the following simple stochastic harmonic resonator 
model (see, e.g., [4, 6]): 



dt2 



i^+^'m-m, (1) 



where £,{t) is temporally white noise, 7 is the damping 
coefficient, and the resonator frequency is defined by the 
angular velocity w (rad/s). Letting the oscillation fre- 
quency change over time and including harmonics allows 
the modeling of more complicated periodic and quasi- 
periodic (almost periodic) properties (c/. [4, 7]). 

However, the oscillatory phenomena can also contain 
spatial properties. This leads to a space-time model, 
where the process can be described by a spatial field 
that is evolving in time. The main contribution herein 
is to set up a model in which the temporal behavior has 
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oscillatory characteristics, such that the process can be 
described by a spatio-temporal resonator model: 



5*2 



dt 



+ Bf{x,t)=^{x,t), (2) 



where A and B are linear operators modeling the space- 
time interactions in the oscillator state. Additionally, 
it is shown how infinite-dimensional Kalman filters and 
smoothers provide computationally effective means of 
computing the Bayesian solution for detecting and es- 
timating the oscillations in noisy measurement data. 

Previously, for discrete space and time, the spatio- 
temporal interactions in such data have been modeled, 
for example, with seasonal VARMA (vector autoregres- 
sive moving average) (see, e.g., [8, 9]) models. However, 
incorporating the spatial structure and predicting new 
measurements in space and time is difficult, if not impos- 
sible, in these models. For continuous-valued treatment, 
one can resort to neural networks (see, e.g., [10]), which 
make it possible to account for the latent structure, but 
provide few tools for assessing the model structure or 
interpreting the results. 

It would also be possible to formulate the spatio- 
temporal model and the related Bayesian estimation 
problem directly in terms of Gaussian processes (GPs), 
for example by using periodic covariance functions [11]. 
Unfortunately, the direct use of this approach leads to an 
intractable cubic computational complexity 0{T^) in the 
number of time steps T. To some extent, it is possible to 
reduce this problem by using sparse approximations (see, 
e.g., [11, 12]), but this does not solve the problem fully. 

The use of stochastic partial differential equation 
(SPDE) based models to form computationally efficient 
solutions to Gaussian process regression problems (or 
equivalent Kriging problems) has recently been discussed 
in [13, 14]. In particular, in [14] the authors propose 
a method for converting a covariance function based 
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spatio-temporal Gaussian process regression model into 
an equivalent SPDE type model. The advantage of this 
approach is that the Bayesian inference problem of the 
resulting model can be solved using infinite-dimensional 
Kalman filtering and smoothing with linear computa- 
tional complexity in time. 

This work follows an approach similar to [14], except 
that the SPDE model is formulated directly as a linear 
combination of spatio-temporal oscillators, rather than 
first forming a covariance function-based Gaussian pro- 
cess regression problem and then converting it into an 
SPDE. Although some modeling freedom is lost in the 
present approach, the advantage is that it always pro- 
duces a model that can be solved with a linear time com- 
plexity algorithm, and no additional conversion proce- 
dures are required to achieve this. 

The following sections introduce the spatio-temporal 
resonator model and explain how to select the spatial 
operators. The operators are chosen such that an or- 
thogonal basis for the model can be formed by using 
the eigenvalue decomposition of the Laplace operator. 
The Hilbert space method approach to solving the GP 
model using Kalman filtering is discussed in brief to- 
gether with maximum likelihood parameter estimation. 
As a proof-of-concept demonstration a one-dimensional 
example is presented. The method is also applied to em- 
pirical weather data (on a spherical surface) and brain 
data (in a polar 2D domain) . 



II. METHODS 

A. Spatio- Temporal Resonator Model 

A model for a general oscillatory phenomenon is con- 
structed as a superposition of several resonators (sepa- 
rate resonators and their harmonics) with known angu- 
lar velocities ujj {i.e. frequencies), but unknown phases 
and amplitudes. These are modeled as spatially inde- 
pendent realizations of stochastic processes. The sum 
^^j^ fj{x,t) of the oscillatory components fj{x,t) can 
be defined using separate state space models. The spa- 
tially independent version of such a resonating field can 
be presented as a partial differential equation (see Eq. (1), 
or [4] for details): 



5<2 



-7r 



dt 



where a; e il (for some domain fl C M") denotes the spa- 
tial variable and t £ represents time. The perturba- 
tion term ^j(x,t) is white noise, both spatially and tem- 
porally . The above formulation also contains a damping 
factor 7j, which was assumed to be zero in [4]. 

Here, this formulation is extended to account for spa- 
tial structure by assuming that the local derivative de- 
pends not only on time, but also on surrounding locations 
through some spatial linear operator. Including linear 



operators that affect both the oscillation and damping 
results in: 



a<2 



This model contains three types of spatial dependency. 
The selection of operators Aj and Bj allows the suit- 
able definition of spatial coupling through the first and 
second temporal derivative. Some spatial and temporal 
structure can also be assumed in the process noise term 
S,j{x,t) through a correlation structure: 

Cj{x,x') ^Y.[£,j{x,t)£,j{x' ,t')] ^C(^^j{x,x')5{t-t'). 



B. Choosing Spatial Operators 

If the operators Aj and Bj are assumed to be trans- 
lation and time invariant, the corresponding Fourier do- 
main transfer functions Aj{iUx) and Bj{iVx) can be cal- 
culated. Taking both spatial and temporal Fourier trans- 
forms of Eq. (3) results in: 

{ivtf'Fj{iu,^,ivt) + {ivt)Aj{iu,^)Fj{iu,^,ivt)+ 



Bj(iVx)Fj(iUx,ivt) = ^jiiVx,m)- 



Solving Fj from above provides: 



Fj{iUx,Wt) 



{ivtY + {wt)Aj{iUx) + Bj{iUj;) ' 
which corresponds to the spectral density: 



Sj{ii>x,m) 



(wt)2 + {wt)Aj{iVx) + Bj{iUx) 



where Qj{vx) = \'^j{'ii^x,ii^t)\'^ is the spectral density of 
£_j. If the operators Aj and Bj are assumed to be for- 
mally Hermitian, the identities Aj{ii>x) = Aj{—iUx) and 
Bj{iVx) — Bj{—Wx) hold, which simplifies the spectral 
density to: 



Sj{iVx,wt) 



Q]{i^x) 



[v^-Bj{iUx)] +iyfAj{iUx) 



The divisor derivative zeros of the system suggest that 
the system has a temporal resonance of I'f = BjiiVx) — 
A^j{wx)/2. The temporal oscillation is included in the 
angular velocity of ojj by setting Bj{iVx) = A|(ii/j,)/2 + 
This gives the spectral density in the form: 



Sj{Ux,Vt) 



According to Bochner's theorem [15] every positive def- 
inite function is the Fourier transform of a Borel measure. 
This requires the spectral density to be positive every- 
where in order to be a valid Fourier transform of a co- 
variance function. This condition is fulfilled MQjivx) is a 
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positive function (i.e. a valid spectral density). To ensure 
the causality and stability of the system Aj^iUx) must be 
chosen such that it is a positive function, which corre- 
sponds to the operator Aj being positive (semi)definite. 
The operator Bj is also chosen to be positive, which re- 
sults in the condition A|(ii>'^)/2 + ujj >0. This holds, if 
Aj is real and positive. Zero values in the spectrum cor- 
respond to infinite peaks. However, this does not seem 
to be a problem, because if both operators are zero the 
model falls back to being spatially independent, where 
the only spatial structure comes from the process noise 
term ^(a;, t). 

To make the model actually useful, some choices must 
be made. The operator Aj must be positive semidefinite. 
Examples of such operators are the identity operator I 
and the negative Laplacian —A = — V'^. Therefore, the 
following operator structure is considered: 



= 2(7, 



where jj,Xj ^ ^-''^ some non- negative constants and 
V* is the so-called biharmonic operator. 

A covariance function for ^j(x,t) must be chosen, 
which can be virtually any spatial stationary covariance 
function 

E[^j{x,t)^j{x',t')] = C^j{x~x')5{t~t'). 

In theory, the covariance function C^j could also be non- 
stationary. 



C. Modeling Spatio- Temporal Data 

Combining all the components in the model provides 
the solution as a superposition of all the oscillator com- 
ponents f{x,t) = X^jLi /j (^i The oscillator compo- 
nent fj{x, t) is defined by a stochastic partial differential 
equation with the Dirichlct boundary conditions: 



+ A. 



dfj{x,t) 
dt 



+ Bjfjix,t) = ^jix,t), 



for (a;, t) e^x M+, and fj{x, t) = for {x, t) edQx M+, 
for all j = 1, 2, . . . , iV. The state of the system is defined 
as a combination of the periodic oscillating fields and 
their first temporal derivatives: 

fix,t) = iMx.t) §^h{x,t) ... fN{x,t) yN[x,t)Y. 

This leads to the linear state space model, which can be 
expressed in the following form: 



df{x,t) 



dt 



Tf{x,t) + L^{x,t) 



(4) 



where is a block-diagonal matrix such that each block 
J consists of a 2 X 2 matrix of linear operators, and L is 
a block-diagonal matrix consisting of 2 x 1 blocks: 



I 

-Bj —Aj 



and L-i 



The measurement model is constructed by defining a 
linear operator fik through which the model is observed 
at discrete time steps t/. at known locations x°^^ G n,i — 
1,2,..., dk. The measurement noise term ^ A/'(0, Rk) 
in Eq. (4) is a Gaussian random variable of dimension 
dk- In step fc, the observed values are yk £ M.'^'' . For 
notational convenience, the possibility of depending 
on time has been omitted. However, it is included in one 
of the demonstrations, where the oscillation frequencies 
in ujj{t) change over time. 



D. Infinite-Dimensional Kalman Filtering 

Eq. (4) is the infinite-dimensional counterpart of a 
continuous-time state space model, where the linear ma- 
trix evolution equation has been replaced by a linear dif- 
ferential operator equation (c/. [14]). The first equation 
(the dynamic model) in (4) is an infinite-dimensional lin- 
ear stochastic differential equation [15]. Here, ^ is a dif- 
ferential operator, and the equation is an evolution type 
stochastic partial differential equation (SPDE) [15, 16]. 

Treating the temporal variable separately in the evo- 
lution type SPDE enables the use of infinite-dimensional 
optimal estimation methods. However, these methods 
are meant for discrete time estimation, and therefore the 
evolution equation needs to be discretized with respect 
to time. First, the evolution operator is formed: 

U{At) = exp(At:F) , 

where exp(-) is the operator exponential function. A so- 
lution to the stochastic equation can now be given as (see 
[14, 15] for details): 

f{x, tk+i) = U{tk+i - tk) f{x, tk) 



Uitk + l-T)L^{x,T)dT, (5) 



where tk+i and tk < ifc+i a-re arbitrary. The 
second term is a Gaussian process with covari- 
ance function Q{x,x';tk,tk+i) = J^^^^ Vl{tk+i — 
t) L C^{x, x') L^U* [tk+i - t) dr. This leads to the fol- 
lowing discrete-time model: 



f{x,tk)^U{Atk)f{x,tk-i) + qk{x) 
Vk ^^kf{x,t) +rk, 



(6) 



where Atk = tk — ife-i and the process noise qk{x) ^ 
Q'P{QtQ{x,x' ; Atk)). This discretization is not an ap- 
proximation, but is the so-called mild solution to the 
infinite-dimensional differential equation [15]. 
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1. Infinite- Dimensional Kalman Filter and Smoother 

The infinite- dimensional Kalman filter [17-19] is a 
closed-form solution to the infinite-dimensional linear fil- 
tering problem (6). Here, a two-step scheme is presented, 
which first calculates the marginal distribution of the 
next step using the known system dynamics. The fol- 
lowing formulation uses a notation similar to [14] and 
can be compared to the finite-dimensional Kalman filter 
[20]. 

The infinite-dimensional prediction step can be ex- 
pressed as follows: 

"^fe|fe-i(a;) =U{Mk)mk-i\k^i{x) 
Ckik-i{x,x')=U{Atk)Ck-iik-i{x,x')U*{Atk) (7) 
+ Q{x,x';Atk), 

where (•)* denotes an adjoint, which in practice swaps the 
roles of inputs x and x' , and operates from the right. The 
operator adjoint can be seen as an operator version of a 
matrix transpose. The recursive iteration is initialized by 
presenting the prior information in the form f{x,to) ~ 
gp {Tnoix),Coix,x')). 

The algorithm then uses each observation to update 
the distribution to match the new information obtained 
by the measurement in step k. This is the infinite- 
dimensional update step: 

Sk = fik Ck\k-i{x, x') Til + Rk 
Kk{x) = Ckik-i{x,x')'HlS-' 
mk\k{x) = mk\k-i{x) + Kk{x) {yk - Tikm^k-iix)) 
Ck\k{x,x') = Ck\k-iix,x') - Kk{x) SkKl{x), 

(8) 

where (■)~^ denotes the matrix inverse. As a result, the 
filtered forward-time posterior process in step k (time tk) 
is given by fk\k{x) ~ GV {mk\k{x), Ck\k{x,x')). 

The purpose of optimal (fixed-interval) smoothing is to 
obtain in closed-form the marginal posterior distribution 
of the state fk in time step tk, which is conditional on 
all the measurements yi-.r, where fc g [1, . . . , T] is a fixed 
interval. 

The infinite-dimensional Rauch-Tung-Striebel (RTS) 
smoother equations are written so that they utilize the 
Kalman filtering results ink\k{x) and Ck\k{x,x') as a 
forward sweep, and then perform a backward sweep to 
update the estimates to match the forthcoming observa- 
tions. The smoother's backward sweep may be expressed 
with the following infinite-dimensional RTS smoothing 
equations [14]: 

rrik+iikix) =U{Atk)mk\k{x) 
Ck+i\k{x, x') = U{Atk) Ck\k{x, x') U*{Atk) 
+ Qk{x,x';Atk) 
Gk = Ck\k{x,x')U*{Atk) [Ck+i\k{x,x')Y' 
mk\T{x) = mk\kix) 



+ Qk [mk+i\T{x) - mk+i\k{x)] 
Ck\T{x, x') = Ck\k{x, x') + Qk {Ck+i\T{x, x') 
- Ck+i\k{x,x'))gl. 

(9) 

In the above equations, (•)^^ denotes the operator in- 
verse. In addition, note that Qk is a linear operator whose 
kernel is defined via the covariance kernels of the filtering 
results. This makes the notation slightly more challeng- 
ing. 

Once both the Kalman filtering and Rauch-Tung- 
Striebel sweeps are performed on the model given the 
observed data, the marginal posterior is obtained, which 
can be represented as the Gaussian process: 

f{x, tk 1 Vi-.t) ^ GV {mk\T{x), Ck\T{x, x')) , 

where the observed values yk G are given at discrete 
time points t^, fc = 1, 2, . . . , T, and measured at known 
locations xf^i^ e fi, i = 1, . . . , d^. The resulting process 
functions can be evaluated at any test point x^, £ Q, 
by simply considering an appropriate measurement func- 
tional T-L. Thus, the marginal posterior of the value of 
f{x^,tk) in cc* at time instant tk is: 

p{fix^,tk) 1 yi-.r) = 

J^{f{x*,tk) 1 mk\T{x*),Ck\T{x^,x^)) . 

Values could also be predicted at more time steps. A 
test time point t^ should be taken into account when 
performing the time discretization. The state of the sys- 
tem f{x,tsf) would be predicted in this step, but because 
there are no data, no updating step is needed. 

One detail worthy of note is the connection between 
the standard (in this case) spatial GP model and the evo- 
lution type state space SPDE. If the temporal evolution 
model is left out, that is ^ = and Qc{x,x') = are 
chosen, the estimation task for a spatial GP model could 
be solved by considering only one measurement step and 
using the same equations. 

2. Hubert Space Methods 

The infinite-dimensional Kalman filtering and smooth- 
ing equations can be converted to a tractable form by 
introducing Hilbert space methods (basis function ap- 
proximations) . The eigenfunction expansion of the linear 
operator is considered and combined with the infinite- 
dimensional framework. By truncating the expansion, 
a finite-dimensional approximate solution is obtained, 
which can be evaluated. 

In the spatio-temporal resonator model, the nega- 
tive Laplace operator can be considered in some spa- 
tial domain f2 subject to Dirichlet boundary conditions. 
This results in an eigenfunction equation in the form 
—S/^ipnix) = Xnipn{x), where ipn{x) is an eigenfunc- 
tion and A„ is the corresponding eigenvalue for each 



5 



n = 1, 2, . . . , A'' and spatial coordinate x G fl. The solu- 
tion presented here f{x,t) will be transformed to a new 
basis, which is given by the eigendecomposition of the 
linear operator in Q. This new basis decodes the spa- 
tial structure so that only f{t) remains, being a finite- 
dimensional approximation of f{x,t). 

After the time-discretization step in (5), the finite- 
dimensional approximation of the s component model 
leads to a sA^-dimensional state space model: 

fk+l = fk + Qk 

Vk = Hkfk + rk, 

where the sparse dynamic model Ak and the noise term 
qk ~ A/'(0, Qk) are given in the basis defined by the eigen- 
function expansion of in (see the supplementary ma- 
terial in [14] for detailed equations). 

E. Parameter Estimation 

All the parameters needed in the model (typically the 
unknowns in the covariance function, the damping con- 
stants, and the measurement noise variance) are sum- 
marized as a vector quantity 6. For notational conve- 
nience, the parameters have not been explicitly written 
out in (7) and (8), but could be included as tl{At,6), 
Q{x,x';t,t'-6) and Rk{0). Based on these results, the 
marginal likelihood of the measurements can be com- 
puted, given 6 (see, e.g., [21]): 

T 

p{yi, ...,yT\e)^l[U{yk\nk mkik-i{x, 6), Sk{e)) . 

k=l 

Hence, the marginal log-likelihood function for maximum 
likelihood (ML) estimation can be given as: 

T T 

m = -\ Elog |2^Sfe(0)|-^ m,|fe_i(aj, 0))^ 

k=l k=l 

X S^\e) {yk-'Hkmk\k-i{x,e)) . (10) 

The parameter estimation problem now reverts to the 
maximization of the log-likelihood function: 6 = 
arg maxe £{9). 

Using the log-likelihood function, the un-normalized 
posterior distribution could also be formed easily, which 
would allow the computation of maximum a posteriori 
(MAP) estimates or the use of a Metropolis-Hastings 
type of Markov chain Monte Carlo (MCMC) for integra- 
tion over the parameters. 

III. RESULTS 
A. Illustrative One-Dimensional Example 

Visualizations for simulated data in one spatial dimen- 
sion are shown in Fig. 1 as an example of the spatio- 
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FIG. 1. This figure illustrates an example of a stochastic res- 
onator realization in one spatial dimension. The left-hand 
figure shows the observation locations and the correspond- 
ing estimate of the oscillating field. The state at t = 0.50 s 
is shown in the right-hand figure, where the estimate is pre- 
sented together with the dashed true values of the process. 

temporal resonator model. The model contains one res- 
onator f{x,t) oscillating at a frequency of 6 Hz, where 
X e [—L,L] over a time-span t G [0, 1] (in seconds). 

The Matern covariance function is considered (see, e.g., 
[11]) for the perturbing dynamic noise. This class of 
stationary isotropic covariance functions is widely used 
in many applications and their parameters have under- 
standable interpretations. A Matern covariance function 
can be expressed as: 

C(.) = .^1^(72^9^.(72^9, (11) 

where r = \\x — x'\\, F(-) is the Gamma function and 
Ki^{-) is the modified Bessel function. The covariance 
function is characterized by three parameters: a smooth- 
ness parameter v, a distance scale parameter I, and a 
strength (magnitude) parameter cr, all of which are pos- 
itive. 

For simulating data, the model parameters were cho- 
sen so that 7 = 1 and x = 0.01. The perturbing dynamic 
noise covariance function parameters were: v = 3/2, 
I = 0.1 L, and s = 25. The Gaussian measurement 
noise variance was cr^ = 0.1^. A truncated eigenfunction 
decomposition with 32 eigenfunctions was used. Alto- 
gether, 2500 noisy observations were considered. 

The model parameters are fitted by optimizing the 
marginal log-likelihood (10) using a conjugate gradient 
optimizer and square root versions of the filtering equa- 
tions (7)-(9) for numerical stability. To avoid bad local 
minima, ten random restarts were attempted, and the 
run with the best marginal likelihood was selected. The 
optimized parameters were ct^ « 0.098^, / sa 0.106L, 
s « 30.8, 7 « 0.690, and x ~ 0.015. 

Fig. 1 shows the estimation mean m(x, t) of the true 
space-time oscillating field f{x,t) as a color surface plot. 
The measurement locations in space-time are shown as 
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(a) Slow temperature bias (b) Oscillatory structure with two harmonics 

FIG. 2. Slow bias and resonator maps for temperatures as a snapshot on July 8, 2011 at 2 PM (GMT). The weather stations 
are marked with crosses and areas of uncertainty are hatched. 



crosses. The oscillatory behavior is clear, and even 
clearer along the slice x = 0, which is shown in the lower 
figure. A slice f{x,0.5) along the spatial dimension is 
shown on the right. Both slices contain the true values 
for / (dashed), the estimate mean (solid), and a shaded 
95% uncertainty interval. 



B. Weather Data on the Surface of the Globe 

In this section, the proposed methods are applied to 
hourly observations of temperature readings in centi- 
grade, which were collected worldwide by National En- 
vironmental Satellite, Data, and Information Services 
(NESDIS)[22]. A subset of the data was considered 
consisting of hourly temperature measurements for one 
month (July, 2011), resulting in a time series of 745 tem- 
poral points. The temperatures were recorded at 11 344 
different spatial locations, for which the longitudinal and 
latitudinal coordinates are known and assumed exact. 
The locations of the stations are marked as crosses in 
Fig. 2a. However, not all stations provide hourly data; 
altogether there are 5 637501 measurements. 

A three-component setup was used, where the first la- 
tent component, a bias term, accounts for the slow drift- 
ing of the mean temperature, and the other two are oscil- 
latory components for the daily variation of temperature. 
The first of these two oscillates at the constant base fre- 
quency of 1/day, and the second is the first harmonic 
(frequency 2/day). The bias term was constructed as an 
oscillator with zero frequency, which can be seen as a 
spatio-temporal Wiener velocity model. The spatial co- 
variance function of the process noise was fixed to the 
squared-exponential covariance function. 

The model parameters (covariance function magni- 
tudes, length scales, and the Gaussian measurement noise 
variance) were all optimized with respect to marginal 
likelihood using a few random restarts to avoid bad lo- 
cal minima. The damping parameters were virtually zero 
in all runs, therefore they were fixed to zero in the final 
model. This means that the spatial dependencies stem 



from the perturbation structure alone. 

The estimation results for the temperature oscillation 
model are presented here as a snapshot of the tempera- 
ture map over the globe on July 8, 2011 at 2 PM (GMT). 
The results in Fig. 2 are split into two in order to show 
the influence of the slow-moving bias and the resonat- 
ing part. Fig. 2a also shows the spatial locations of the 
weather stations. The hatched regions indicate uncer- 
tainty (standard deviation > 2°C), which corresponds 
to the regions with very few or no observations. 

This test setup is subject to many simplifications and 
assumptions, which affect the results; the surface of the 
Earth is actually not a symmetrical sphere and the evi- 
dent fact that the fluctuation covariance structure is not 
stationary is disregarded. However, the setup clearly 
captures two effects: the summer on the northern hemi- 
sphere and the day-night variation (afternoon in Europe 
and Africa in the figure). 



C. Modeling Pulsations in the Brain 

Recent advances in functional magnetic resonance 
imaging (fMRI) techniques have demonstrated the im- 
portance of computational methods in modeling brain 
data. In [4], it was shown that eliminating oscillating 
physiological noise components in fMRI can be achieved 
by using a spatially independent resonator model com- 
bined with Kalman filtering. This method was named 
DRIFTER. This approach is now extended by showing 
how to account also for spatial dependencies. 

One ~ 30 s run of empirical fMRI data was considered. 
The set of fast-sampled data of one slice is used here to 
demonstrate the spatio-temporal resonator model in two- 
dimensional polar coordinates. This fMRI data, together 
with anatomical images for one volunteer, were obtained 
using a 3T scanner (Siemens Skyra) located at the Ad- 
vanced Magnetic Imaging Centre (AMI) of Aalto Uni- 
versity School of Science using a 32-channel receive-only 
head coil. For the functional imaging, the major param- 
eters were repetition time (TR) 77 ms, echo time (TE) 
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(a) Cardiac noise amplitude (b) Respiratory noise amplitude 

FIG. 3. Mean amplitude maps for the physiological noise components. The results are shown both in isolation and overlaid on 
top of the corresponding anatomical image. 



20 ms, flip angle (FA) 60°, field-of-view (FOV) 224 mm, 
and matrix size 64x64. The measurements were per- 
formed as part of AMI Centre's local technical methods 
development research and conformed to the guidelines 
of the Declaration of Helsinki. The research was ap- 
proved by the ethical committee in the Hospital District 
of Helsinki and Uusimaa. 

In order to simulate resting state conditions, the stim- 
ulus was a fixed dot in the center of the visual field of 
the volunteer. The heart and respiratory signals were 
recorded, and time-locked to the fMRI data during the 
run. The oscillation frequency of the physiological noise 
components was quasi-periodic (rather than exactly pe- 
riodic), which implies that the frequencies change over 
time. External reference signals with the interacting mul- 
tiple model (IMM) approach, presented in [4], were used 
to estimate the frequency time series of heart beats and 
respiration cycles. The cardiac frequency alternated be- 
tween 54-64 beats per minute and the respiratory fre- 
quency between 6-15 cycles per minute. 

One slice of fast-sampled fMRI data was used. The 
sampling interval was 0.077 s and the whole 64x 64 matrix 
was observed during each of 350 time steps. The spatial 
domain Q, was chosen to be a circular disk with a radius 
of « 155 mm. The spatio-temporal resonator model has 
three components: a slowly moving brain blood-oxygen- 
level-dependent (BOLD) signal — which also includes 
scanner drift and other slow phenomena — modeled with 
a spatio-temporal Wiener velocity model (see Sec. IIIB), 
and two space-time resonators oscillating at the time- 
dependent cardiac and respiratory frequencies, respec- 
tively. Here, only the resonators for the base frequencies 
were included; more complex signals could be accounted 
for by including harmonics. An eigenfunction decompo- 
sition of the linear operators in f2 with 300 eigenfunctions 
was used for each of the three components. Model param- 
eters were chosen by studying the spatially independent 
model first. 

Figure 3 shows the spatial amplitude of physiological 
noise contribution averaged over time. The results re- 
semble the maps in [5], where the oscillators were as- 
sumed to be spatially independent and only the final re- 



sults were spatially smoothed. This suggests that the 
method is able to capture the space-time structure of 
the oscillations. The cardiac influence is strong near the 
large cerebral arteries (see Fig. 3a), and the respiration 
causes artifacts near the eyes that are partly induced by 
movement. 



IV. CONCLUSION AND DISCUSSION 

This paper proposed a computationally effective 
stochastic partial differential equation model and an in- 
ference method for detecting and modeling latent os- 
cillatory structures in spatio-temporal data. It showed 
how a physical first-principles SPDE model for spatio- 
temporal oscillations can be constructed, and how 
the Bayesian inference can be effectively applied using 
infinite-dimensional Kalman filtering and Hilbert space 
methods. This filtering is related to Gaussian process 
regression and Hilbert space valued stochastic processes. 
The proposed method allows a reduction of the complex- 
ity of a direct GP solution from cubic to linear with re- 
spect to measurements in the temporal dimension. 

A truncated eigenfunction expansion of the Laplace 
operator was used to form a finite-dimensional basis over 
the spatial domain, which made it possible to revert to 
the traditional Kalman filtering scheme. The eigenfunc- 
tion expansions of the Laplace operator in both spherical 
and Cartesian coordinates were used in the numerical 
computations. 

The numerical results show that the truncated ex- 
pansion puts some restrictions on the spatial short-scale 
variability. The basis function approach tends to make 
the model spatially smooth, a problem that has been 
dealt with before in many ways under the GP regression 
scheme (see, e.(?., [11]). However, in many applications 
such as in functional brain data analysis, this is not a 
problem since a few hundred basis functions are suffi- 
cient to match the required spatial resolution. Several 
methodological extensions could be considered as well, 
for example the possibility of including non-stationary 
covariance functions in the process noise model. 
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In general, the spatio-temporal model can mitigate the 
problems related to slow sampling rates, because the spa- 
tial information can compensate for missing temporal 
data. This turns such models into powerful tools for sig- 
nal reconstruction and noise elimination, for example in 
fMRI studies, especially as the computational complexity 
grows only linearly with the length of the measurement 
session. 
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